<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>cggbal.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="CGGBAL.1"></a><a href="cggbal.f.html#CGGBAL.1">CGGBAL</a>( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
     $                   RSCALE, WORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      CHARACTER          JOB
      INTEGER            IHI, ILO, INFO, LDA, LDB, N
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      REAL               LSCALE( * ), RSCALE( * ), WORK( * )
      COMPLEX            A( LDA, * ), B( LDB, * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="CGGBAL.20"></a><a href="cggbal.f.html#CGGBAL.1">CGGBAL</a> balances a pair of general complex matrices (A,B).  This
</span><span class="comment">*</span><span class="comment">  involves, first, permuting A and B by similarity transformations to
</span><span class="comment">*</span><span class="comment">  isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N
</span><span class="comment">*</span><span class="comment">  elements on the diagonal; and second, applying a diagonal similarity
</span><span class="comment">*</span><span class="comment">  transformation to rows and columns ILO to IHI to make the rows
</span><span class="comment">*</span><span class="comment">  and columns as close in norm as possible. Both steps are optional.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Balancing may reduce the 1-norm of the matrices, and improve the
</span><span class="comment">*</span><span class="comment">  accuracy of the computed eigenvalues and/or eigenvectors in the
</span><span class="comment">*</span><span class="comment">  generalized eigenvalue problem A*x = lambda*B*x.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  JOB     (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">          Specifies the operations to be performed on A and B:
</span><span class="comment">*</span><span class="comment">          = 'N':  none:  simply set ILO = 1, IHI = N, LSCALE(I) = 1.0
</span><span class="comment">*</span><span class="comment">                  and RSCALE(I) = 1.0 for i=1,...,N;
</span><span class="comment">*</span><span class="comment">          = 'P':  permute only;
</span><span class="comment">*</span><span class="comment">          = 'S':  scale only;
</span><span class="comment">*</span><span class="comment">          = 'B':  both permute and scale.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The order of the matrices A and B.  N &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  A       (input/output) COMPLEX array, dimension (LDA,N)
</span><span class="comment">*</span><span class="comment">          On entry, the input matrix A.
</span><span class="comment">*</span><span class="comment">          On exit, A is overwritten by the balanced matrix.
</span><span class="comment">*</span><span class="comment">          If JOB = 'N', A is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDA     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array A. LDA &gt;= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  B       (input/output) COMPLEX array, dimension (LDB,N)
</span><span class="comment">*</span><span class="comment">          On entry, the input matrix B.
</span><span class="comment">*</span><span class="comment">          On exit, B is overwritten by the balanced matrix.
</span><span class="comment">*</span><span class="comment">          If JOB = 'N', B is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDB     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array B. LDB &gt;= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  ILO     (output) INTEGER
</span><span class="comment">*</span><span class="comment">  IHI     (output) INTEGER
</span><span class="comment">*</span><span class="comment">          ILO and IHI are set to integers such that on exit
</span><span class="comment">*</span><span class="comment">          A(i,j) = 0 and B(i,j) = 0 if i &gt; j and
</span><span class="comment">*</span><span class="comment">          j = 1,...,ILO-1 or i = IHI+1,...,N.
</span><span class="comment">*</span><span class="comment">          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LSCALE  (output) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">          Details of the permutations and scaling factors applied
</span><span class="comment">*</span><span class="comment">          to the left side of A and B.  If P(j) is the index of the
</span><span class="comment">*</span><span class="comment">          row interchanged with row j, and D(j) is the scaling factor
</span><span class="comment">*</span><span class="comment">          applied to row j, then
</span><span class="comment">*</span><span class="comment">            LSCALE(j) = P(j)    for J = 1,...,ILO-1
</span><span class="comment">*</span><span class="comment">                      = D(j)    for J = ILO,...,IHI
</span><span class="comment">*</span><span class="comment">                      = P(j)    for J = IHI+1,...,N.
</span><span class="comment">*</span><span class="comment">          The order in which the interchanges are made is N to IHI+1,
</span><span class="comment">*</span><span class="comment">          then 1 to ILO-1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RSCALE  (output) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">          Details of the permutations and scaling factors applied
</span><span class="comment">*</span><span class="comment">          to the right side of A and B.  If P(j) is the index of the
</span><span class="comment">*</span><span class="comment">          column interchanged with column j, and D(j) is the scaling
</span><span class="comment">*</span><span class="comment">          factor applied to column j, then
</span><span class="comment">*</span><span class="comment">            RSCALE(j) = P(j)    for J = 1,...,ILO-1
</span><span class="comment">*</span><span class="comment">                      = D(j)    for J = ILO,...,IHI
</span><span class="comment">*</span><span class="comment">                      = P(j)    for J = IHI+1,...,N.
</span><span class="comment">*</span><span class="comment">          The order in which the interchanges are made is N to IHI+1,
</span><span class="comment">*</span><span class="comment">          then 1 to ILO-1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace) REAL array, dimension (lwork)
</span><span class="comment">*</span><span class="comment">          lwork must be at least max(1,6*N) when JOB = 'S' or 'B', and
</span><span class="comment">*</span><span class="comment">          at least 1 when JOB = 'N' or 'P'.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0:  successful exit
</span><span class="comment">*</span><span class="comment">          &lt; 0:  if INFO = -i, the i-th argument had an illegal value.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  See R.C. WARD, Balancing the generalized eigenvalue problem,
</span><span class="comment">*</span><span class="comment">                 SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      REAL               ZERO, HALF, ONE
      PARAMETER          ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0 )
      REAL               THREE, SCLFAC
      PARAMETER          ( THREE = 3.0E+0, SCLFAC = 1.0E+1 )
      COMPLEX            CZERO
      PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ) )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      INTEGER            I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1,
     $                   K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN,
     $                   M, NR, NRP2
      REAL               ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2,
     $                   COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX,
     $                   SFMIN, SUM, T, TA, TB, TC
      COMPLEX            CDUM
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.124"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      INTEGER            ICAMAX
      REAL               SDOT, <a name="SLAMCH.126"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
      EXTERNAL           <a name="LSAME.127"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, ICAMAX, SDOT, <a name="SLAMCH.127"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           CSSCAL, CSWAP, SAXPY, SSCAL, <a name="XERBLA.130"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, AIMAG, INT, LOG10, MAX, MIN, REAL, SIGN
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Functions ..
</span>      REAL               CABS1
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Function definitions ..
</span>      CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Test the input parameters
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      IF( .NOT.<a name="LSAME.146"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'N'</span> ) .AND. .NOT.<a name="LSAME.146"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'P'</span> ) .AND.
     $    .NOT.<a name="LSAME.147"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'S'</span> ) .AND. .NOT.<a name="LSAME.147"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'B'</span> ) ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -4
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -6
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.157"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CGGBAL.157"></a><a href="cggbal.f.html#CGGBAL.1">CGGBAL</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.0 ) THEN
         ILO = 1
         IHI = N
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.1 ) THEN
         ILO = 1
         IHI = N
         LSCALE( 1 ) = ONE
         RSCALE( 1 ) = ONE
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( <a name="LSAME.177"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'N'</span> ) ) THEN
         ILO = 1
         IHI = N
         DO 10 I = 1, N
            LSCALE( I ) = ONE
            RSCALE( I ) = ONE
   10    CONTINUE
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      K = 1
      L = N
      IF( <a name="LSAME.189"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'S'</span> ) )
     $   GO TO 190
<span class="comment">*</span><span class="comment">
</span>      GO TO 30
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Permute the matrices A and B to isolate the eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Find row with one nonzero in columns 1 through L
</span><span class="comment">*</span><span class="comment">
</span>   20 CONTINUE
      L = LM1
      IF( L.NE.1 )
     $   GO TO 30
<span class="comment">*</span><span class="comment">
</span>      RSCALE( 1 ) = ONE
      LSCALE( 1 ) = ONE
      GO TO 190
<span class="comment">*</span><span class="comment">
</span>   30 CONTINUE
      LM1 = L - 1
      DO 80 I = L, 1, -1
         DO 40 J = 1, LM1
            JP1 = J + 1
            IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
     $         GO TO 50
   40    CONTINUE
         J = L
         GO TO 70
<span class="comment">*</span><span class="comment">
</span>   50    CONTINUE
         DO 60 J = JP1, L
            IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
     $         GO TO 80
   60    CONTINUE
         J = JP1 - 1
<span class="comment">*</span><span class="comment">
</span>   70    CONTINUE
         M = L
         IFLOW = 1
         GO TO 160
   80 CONTINUE
      GO TO 100
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Find column with one nonzero in rows K through N
</span><span class="comment">*</span><span class="comment">
</span>   90 CONTINUE
      K = K + 1
<span class="comment">*</span><span class="comment">
</span>  100 CONTINUE
      DO 150 J = K, L
         DO 110 I = K, LM1
            IP1 = I + 1
            IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
     $         GO TO 120
  110    CONTINUE
         I = L
         GO TO 140
  120    CONTINUE
         DO 130 I = IP1, L
            IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
     $         GO TO 150
  130    CONTINUE
         I = IP1 - 1
  140    CONTINUE
         M = K
         IFLOW = 2
         GO TO 160
  150 CONTINUE
      GO TO 190
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Permute rows M and I
</span><span class="comment">*</span><span class="comment">
</span>  160 CONTINUE
      LSCALE( M ) = I
      IF( I.EQ.M )
     $   GO TO 170
      CALL CSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA )
      CALL CSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Permute columns M and J
</span><span class="comment">*</span><span class="comment">
</span>  170 CONTINUE
      RSCALE( M ) = J
      IF( J.EQ.M )
     $   GO TO 180
      CALL CSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
      CALL CSWAP( L, B( 1, J ), 1, B( 1, M ), 1 )
<span class="comment">*</span><span class="comment">
</span>  180 CONTINUE
      GO TO ( 20, 90 )IFLOW
<span class="comment">*</span><span class="comment">
</span>  190 CONTINUE
      ILO = K
      IHI = L
<span class="comment">*</span><span class="comment">
</span>      IF( <a name="LSAME.284"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'P'</span> ) ) THEN
         DO 195 I = ILO, IHI
            LSCALE( I ) = ONE
            RSCALE( I ) = ONE
  195    CONTINUE
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( ILO.EQ.IHI )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Balance the submatrix in rows ILO to IHI.
</span><span class="comment">*</span><span class="comment">
</span>      NR = IHI - ILO + 1
      DO 200 I = ILO, IHI
         RSCALE( I ) = ZERO
         LSCALE( I ) = ZERO
<span class="comment">*</span><span class="comment">
</span>         WORK( I ) = ZERO
         WORK( I+N ) = ZERO
         WORK( I+2*N ) = ZERO
         WORK( I+3*N ) = ZERO
         WORK( I+4*N ) = ZERO
         WORK( I+5*N ) = ZERO
  200 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute right side vector in resulting linear equations
</span><span class="comment">*</span><span class="comment">
</span>      BASL = LOG10( SCLFAC )
      DO 240 I = ILO, IHI
         DO 230 J = ILO, IHI
            IF( A( I, J ).EQ.CZERO ) THEN
               TA = ZERO
               GO TO 210
            END IF
            TA = LOG10( CABS1( A( I, J ) ) ) / BASL
<span class="comment">*</span><span class="comment">
</span>  210       CONTINUE
            IF( B( I, J ).EQ.CZERO ) THEN
               TB = ZERO
               GO TO 220
            END IF
            TB = LOG10( CABS1( B( I, J ) ) ) / BASL
<span class="comment">*</span><span class="comment">
</span>  220       CONTINUE
            WORK( I+4*N ) = WORK( I+4*N ) - TA - TB
            WORK( J+5*N ) = WORK( J+5*N ) - TA - TB
  230    CONTINUE
  240 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      COEF = ONE / REAL( 2*NR )
      COEF2 = COEF*COEF
      COEF5 = HALF*COEF2
      NRP2 = NR + 2
      BETA = ZERO
      IT = 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Start generalized conjugate gradient iteration
</span><span class="comment">*</span><span class="comment">
</span>  250 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      GAMMA = SDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) +
     $        SDOT( NR, WORK( ILO+5*N ), 1, WORK( ILO+5*N ), 1 )
<span class="comment">*</span><span class="comment">
</span>      EW = ZERO
      EWC = ZERO
      DO 260 I = ILO, IHI
         EW = EW + WORK( I+4*N )
         EWC = EWC + WORK( I+5*N )
  260 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2
      IF( GAMMA.EQ.ZERO )
     $   GO TO 350
      IF( IT.NE.1 )
     $   BETA = GAMMA / PGAMMA
      T = COEF5*( EWC-THREE*EW )
      TC = COEF5*( EW-THREE*EWC )
<span class="comment">*</span><span class="comment">
</span>      CALL SSCAL( NR, BETA, WORK( ILO ), 1 )
      CALL SSCAL( NR, BETA, WORK( ILO+N ), 1 )
<span class="comment">*</span><span class="comment">
</span>      CALL SAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 )
      CALL SAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 )
<span class="comment">*</span><span class="comment">
</span>      DO 270 I = ILO, IHI
         WORK( I ) = WORK( I ) + TC
         WORK( I+N ) = WORK( I+N ) + T
  270 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Apply matrix to vector
</span><span class="comment">*</span><span class="comment">
</span>      DO 300 I = ILO, IHI
         KOUNT = 0
         SUM = ZERO
         DO 290 J = ILO, IHI
            IF( A( I, J ).EQ.CZERO )
     $         GO TO 280
            KOUNT = KOUNT + 1
            SUM = SUM + WORK( J )
  280       CONTINUE
            IF( B( I, J ).EQ.CZERO )
     $         GO TO 290
            KOUNT = KOUNT + 1
            SUM = SUM + WORK( J )
  290    CONTINUE
         WORK( I+2*N ) = REAL( KOUNT )*WORK( I+N ) + SUM
  300 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      DO 330 J = ILO, IHI
         KOUNT = 0
         SUM = ZERO
         DO 320 I = ILO, IHI
            IF( A( I, J ).EQ.CZERO )
     $         GO TO 310
            KOUNT = KOUNT + 1
            SUM = SUM + WORK( I+N )
  310       CONTINUE
            IF( B( I, J ).EQ.CZERO )
     $         GO TO 320
            KOUNT = KOUNT + 1
            SUM = SUM + WORK( I+N )
  320    CONTINUE
         WORK( J+3*N ) = REAL( KOUNT )*WORK( J ) + SUM
  330 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      SUM = SDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) +
     $      SDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 )
      ALPHA = GAMMA / SUM
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Determine correction to current iteration
</span><span class="comment">*</span><span class="comment">
</span>      CMAX = ZERO
      DO 340 I = ILO, IHI
         COR = ALPHA*WORK( I+N )
         IF( ABS( COR ).GT.CMAX )
     $      CMAX = ABS( COR )
         LSCALE( I ) = LSCALE( I ) + COR
         COR = ALPHA*WORK( I )
         IF( ABS( COR ).GT.CMAX )
     $      CMAX = ABS( COR )
         RSCALE( I ) = RSCALE( I ) + COR
  340 CONTINUE
      IF( CMAX.LT.HALF )
     $   GO TO 350
<span class="comment">*</span><span class="comment">
</span>      CALL SAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 )
      CALL SAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 )
<span class="comment">*</span><span class="comment">
</span>      PGAMMA = GAMMA
      IT = IT + 1
      IF( IT.LE.NRP2 )
     $   GO TO 250
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End generalized conjugate gradient iteration
</span><span class="comment">*</span><span class="comment">
</span>  350 CONTINUE
      SFMIN = <a name="SLAMCH.441"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'S'</span> )
      SFMAX = ONE / SFMIN
      LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE )
      LSFMAX = INT( LOG10( SFMAX ) / BASL )
      DO 360 I = ILO, IHI
         IRAB = ICAMAX( N-ILO+1, A( I, ILO ), LDA )
         RAB = ABS( A( I, IRAB+ILO-1 ) )
         IRAB = ICAMAX( N-ILO+1, B( I, ILO ), LDB )
         RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) )
         LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE )
         IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) )
         IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB )
         LSCALE( I ) = SCLFAC**IR
         ICAB = ICAMAX( IHI, A( 1, I ), 1 )
         CAB = ABS( A( ICAB, I ) )
         ICAB = ICAMAX( IHI, B( 1, I ), 1 )
         CAB = MAX( CAB, ABS( B( ICAB, I ) ) )
         LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE )
         JC = RSCALE( I ) + SIGN( HALF, RSCALE( I ) )
         JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB )
         RSCALE( I ) = SCLFAC**JC
  360 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Row scaling of matrices A and B
</span><span class="comment">*</span><span class="comment">
</span>      DO 370 I = ILO, IHI
         CALL CSSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA )
         CALL CSSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB )
  370 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Column scaling of matrices A and B
</span><span class="comment">*</span><span class="comment">
</span>      DO 380 J = ILO, IHI
         CALL CSSCAL( IHI, RSCALE( J ), A( 1, J ), 1 )
         CALL CSSCAL( IHI, RSCALE( J ), B( 1, J ), 1 )
  380 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="CGGBAL.480"></a><a href="cggbal.f.html#CGGBAL.1">CGGBAL</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
